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Abstract 

We carry out a completely first-principles study of the ferroelectric phase 
transitions in BaTi 03 . Our approach takes advantage of two features of these 
transitions: the structural changes are small, and only low-energy distortions 
are important. Based on these observations, we make systematically improv¬ 
able approximations which enable the parameterization of the complicated 
energy surface. The parameters are determined from first-principles total- 
energy calculations using ultra-soft pseudopotentials and a preconditioned 
conjugate-gradient scheme. The resulting effective Hamiltonian is then solved 
by Monte Carlo simulation. The calculated phase sequence, transition tem¬ 
peratures, latent heats, and spontaneous polarizations are all in good agree¬ 
ment with experiment. We find the transitions to be intermediate between 
order-disorder and displacive character. We find all three phase transitions 
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to be of first order. The roles of different interactions are discussed. 
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I. INTRODUCTION 


Because of their simple crystal structure, the perovskite oxides present a special oppor¬ 
tunity for the development of a detailed theoretical understanding of the ferroelectric phase 
transition. Within this family of materials, one finds transitions to a wide variety of low- 
symmetry phases, including ferroelectric and antiferroelectric transitions. Both first- and 
second-order transitions are observed, with a full spectrum of transition behavior ranging 
from displacive to order-disorder. The properties of BaTiC> 3 , a much-studied prototypical 
example of this class of compounds, 0 exemplify this rich behavior. BaTiC >3 undergoes a 
succession of first-order phase transitions, from the high-temperature high-symmetry cu¬ 
bic perovskite phase (Fig. f|) to slightly distorted ferroelectric structures with tetragonal, 
orthorhombic and rhombohedral symmetry. There is increasing evidence that the cubic-to- 
tetragonal transition, at first thought to be of the simple displacive kind, may instead be 
better described as of the order-disorder type. 

The variety exhibited by the perovskite oxides shows that the phase transformation 
behavior depends on details of the chemistry and structural energetics of each particular 
compound. Therefore, it is of the first importance to develop a microscopic theory of the 
materials properties which determine the ordering of the phases, the character and thermody¬ 
namic order of the transitions, and the transition temperatures. The value of a microscopic 
approach has long been appreciated, but its realization was hindered by the difficulty of 
determining microscopic parameters for individual compounds. The forms of phenomeno¬ 
logical model HainiltonianJlHi were limited by the available experimental data, leading to 
oversimplification and ambiguities in interpretation. For the perovskite oxides, empirical! 
and non-empirical pair potential methods! did not offer the high accuracy needed for the 
construction of realistic models. 

First-principles density functional calculations offer an attractive approach for enhanc¬ 
ing our microscopic understanding of perovskites and other ferroclectrics. The all-electron 
Full-potential Linearized-Augmented-Plane-Wave (FLAPW) method has been used by sev- 
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eral groups to study ferroelectricity in perovskites within the Local Density Approximation 
(LDA).0’i Recently, King-Smith and Vanderbilt performed a systematic study of structural 
and dynamical properties and energy surfaces for eight common perovskites, ^0 using the 
first-principles ultra-soft pseudopotential method and the LDA. These calculations demon¬ 
strate that ferroelectricity in the perovskites reflects a delicate balance between long-range 
electrostatic forces which favor the ferroelectric state, and short-range repulsions which fa¬ 
vor the cubic phase. While constrained to calculations of zero-temperature properties, these 
calculations yield correct predictions of ground state structures and occurrence of ferroelec¬ 
tric phases for certain materials. They show that high-quality LDA calculations can provide 
considerable insight into the nature of the total-energy surface in the perovskites. For further 
insight into the energetics of ferroelectric compounds, the polarization generated by various 
distortions can be studied directly, using a recent first-principles method by King-Smith 
and Vanderbilt .0 This approach has been applied to the investigation of the zone-center 
phonons in the common perovskite oxides.0 

The application of these first-principles methods can clearly form a foundation for the 
realistic study of the finite-temperature phase transitions. While an ab-initio molecular- 
dynamics simulation of the structural phase transition is not computationally feasible at 
present, we purse an alternative first-principles approach to study ferroelectric phase transi¬ 
tions and demonstrate its application to BaTiC^. In particular, we (i) construct an effective 
Hamiltonian to describe the important degrees of freedom of the system,00 (ii) deter¬ 
mine all the parameters of this effective Hamiltonian from high-accuracy ab-initio LDA 
calculations, 100 and (iii) carry out Monte Carlo (MC) simulations to determine the phase 
transformation behavior of the resulting system. An abbreviated presentation of this work 
has already appeared in Ref. |1^. 

The remainder of this paper is organized as follows. In Sec. II, we go through the detailed 
procedure for the construction of the effective Hamiltonian, and give the explicit formula. In 
Sec. Ill, we describe our first-principles calculations and the determination of the expansion 
parameters in the Hamiltonian. The technical details of the Monte Carlo simulation are 
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presented in Sec. IV. In Sec. V, we report our calculated transition temperatures, order 
parameters, and phase diagram, as well as thermodynamic order and nature of the phase 
transitions. The role of different interactions in determining the phase transition behavior 
is also discussed. Sec. VI concludes the paper. 

II. CONSTRUCTION OF THE HAMILTONIAN 

A. Approximations and local modes 

The central quantity for studying the equilibrium properties of a system at finite tem¬ 
perature is its partition function. This can be determined from the energy surface, i.e., the 
total potential energy as a functional of the atomic coordinates. Since the contribution to 
the partition function decreases exponentially with increasing energy, it is possible to obtain 
an accurate partition function for low-temperature applications from a simplified energy 
surface including only low-energy configurations. Our goal is to construct a parameterized 
Hamiltonian which (i) is ab initio , involving no empirical or semi-empirical input; (ii) results 
in an accurate partition function for the temperature range of interest; (iii) is fully specified 
by a few ab-initio total energy calculations; and (iv) involves only approximations that are 
systematically improvable and removable. 

Our first fundamental approximation is to use an energy surface represented by a low- 
order Taylor expansion. Both experiments and first-principles total energy calculations 
suggest that the ferroelectric (FE) phase transition involves only very small atomic dis¬ 
placements and strain deformations from the equilibrium cubic structure. It is reasonable 
to assume that all the atomic configurations with significant contribution to the partition 
function would be close to this cubic structure in the temperature range of interest. Thus, 
it is natural to represent the energy surface by a Taylor series in the displacements from 
the cubic structure. We include up to fourth-order terms in our expansion; this is clearly 
a minimum, since ferroelectricity is intrinsically an anharmonic phenomenon. By including 
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higher-order terms, this approximation could later be systematically improved. 

It is convenient to describe the small distortions from the cubic structure in terms of the 
three acoustic and twelve optical normal-mode coordinates per k-point. While this could 
be regarded as only a change of basis, it motivates our second fundamental approximation, 
which is to restrict the expansion to include only low-energy distortions. To achieve this sep¬ 
aration, we note that both experimentally-measured and LDA-calculated phonon dispersion 
relations0Ei show that only the lowest TO modes (soft modes) and long-wavelength acoustic 
phonons (strain variables) make significant contributions to the phonon density of states at 
low energy. Experimental studies also suggest that the FE phase transitions are accompa¬ 
nied by a softening of the lowest TO mode and the appearance of a strain. All other phonons 
are hardly affected by the transitions. It is then our second approximation to express the en¬ 
ergy surface only as a function of the soft-mode amplitudes and strain. This approximation 
reduces the number of degrees of freedom per cell from fifteen to six, and greatly reduces 
the number of interaction parameters needed. If necessary, this approximation could later 
be relaxed by including additional modes. 

It is convenient to describe the soft mode over the whole Brillouin zone (BZ) in terms 
of a collective motion of “local modes,” just as one describes an acoustic phonon in terms 
of a collective displacement of individual atoms. However, there is more than one choice of 
local mode which will generate the same soft mode throughout the BZ; an intelligent choice 
simplify the Hamiltonian and reduce the number of calculations needed.0 First, the 
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local mode should be as symmetric as possible, so as to minimize the number of expansion 
parameters needed. Second, the interactions between local modes at different sites are more 
difficult to treat than their on-site energy, so the local mode should be chosen so as to 
minimize intersite interactions. For perovskite ABO 3 compounds, the highest symmetry is 
achieved by centering the local mode on either atom A or B. In the case of BaTi0 3 , the Ti-0 
bond is much stronger than the Ba-0 bond and the motion of the Ti is more important in 
the FE transition, so we choose the local mode which is centered on the Ti atom. 

The soft zone-center (k=0) FE mode in BaTiO.j is a T^ mode which can be characterized 


6 


by the four parameters £4, £#, £o||, and £o_l (for a mode polarized along the j’th Cartesian 
direction, these refer to the displacements of the A atom, the B atom, the O atom that 
forms a B-0 bond along direction j, and the other two O atoms, respectively). We take the 
corresponding local mode to consist of a motion of the central A atom by amount £ 4 , the 
eight neighboring B atoms by amounts £_b/8, and the six neighboring O atoms by amounts 
£o||/2 or £o_l/2, along the j’th Cartesian direction. This mode is illustrated in Fig. flU); its 
amplitude is denoted u r An arbitrary k — 0 soft mode can then be realized as a linear 
superposition of these local modes having identical amplitudes (u x ,u y ,u z ) = u in every cell. 

The harmonic interactions between the local-mode amplitudes u* connecting neighbor¬ 
ing cells i must be chosen to reproduce the harmonic behavior of the soft-mode branch 
throughout the BZ. Long-range Coulomb forces are known to play an important part in 
these interactions; they are characterized in terms of the calculated Born (or “transverse”) 
effective charges.B Thus, the harmonic intersite interactions are represented by a sum of 
two contributions: an infinite-range piece that is precisely the interaction of point dipoles 
whose magnitude is given by the Born effective charge, and corrections which we take to be 
of covalent origin and therefore local. 

To be completely general, anharmonic intercell interactions between neighboring Uj would 
likewise have to be included. Instead, we include only on-site anharmonic interactions, which 
are chosen in such a way that the anharmonic couplings for k — 0 modes of the real system 
are correctly reproduced. This “local anharmonicity approximation” is an important feature 
which helps make our scheme tractable and efficient. To go beyond this approximation, one 
could carry out a careful series of frozen-phonon LDA calculations on supercells to determine 
anharmonic couplings at other points in the BZ. However, past experience has shown that 
calculations of this kind are very cumbersome because of the large number of parameters 
which have to be determined.^ 

With these approximations, our Hamiltonian consists of five parts: a local-mode self 
energy, a long-range dipole-dipole interaction, a short-range interaction between soft modes, 
an elastic energy, and an interaction between the local modes and local strain. Symbolically, 
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E tot = £sdf({ u }) + ^Pl({u}) + £ short ({u}) 
+ £ elas (M) + £ int ({u}, {rn}) , 


( 1 ) 


where u is the local soft-mode amplitude vector, and rft is the six-component local strain 
tensor in Voigt notation (r/i = e n , >H — 2 e 2 3 ). In the following subsections, we present the 
explicit formulae for these five contributions. 


B. Local mode self energy 


The first term is 


£" if ({u})=i:£(u.), <2) 

i 

where E(ui) is the energy of an isolated local mode at cell R, with amplitude u,, relative 
to that of the perfect cubic structure. To describe the FE phase, E(Ui) must contain 
anharmonic as well as harmonic contributions. Since the reference structure is cubic, only 
even-order terms can enter; we choose to truncate at fourth order. Symmetry considerations 
then require that it take the form 

E(ui) = k 2 u 2 + auf + 7 (u 2 ix u 2 iy + u 2 y u 2 z + u 2 z u 2 x ) (3) 

where iq = u,; |, and k 2i a, and 7 are expansion parameters to be determined from first- 
principles calculations. 


C. Dipole-dipole interaction 

The second term in the effective Hamiltonian represents long-range interactions between 
local modes. Only dipole-dipole interactions are considered, since higher-order terms tend to 
be of short range and their effect will be included in the short-range contribution E short ({u}). 
The dipole moment associated with the local mode in cell i is d,; = Z* U;. Here, Z* is the 
Born effective charge for the soft mode, which can be obtained as 



Z* — £,aZ* a + £bZ* b + io\\Z*o\\ + 2£o±Zq ± 


( 4 ) 


from the eigenvector of the soft mode, once the Born effective charges for the ions are 


known. |12| 

In atomic units (energy in Hartree), we have 


Z ~ \ Uj ■ u j 3(■ Uj)(R,jj • uj) 



(5) 


Here, e^, is the optical dielectric constant of the material, Rjj = | R/y |, Rq- = — R,-, and 



In practice, Eq. (|5|) is not directly useful for three-dimensional systems with periodic 
boundary conditions; instead, we use an Ewald construction to evaluate E dpl . We effectively 
terminate the sum in such a way that the k = 0 modes of the supercell will represent physical 
TO(r) modes. For a TO mode, the induced depolarization electric field is zero; from the 
point of view of the dipole sum, it is as though the material were surrounded by a layer of 
metal. In the Ewald construction, this is equivalent to setting the surface terms to zeroJi! 
Under these conditions, and choosing the decay A of the Gaussian charge packets to be small 
enough so that the real-space summation can be entirely neglected, we have 



where is the cell volume and G is the reciprocal lattice vector. 

Because of its long-range nature, the calculation of E dpl is the most time-consuming 
part of our Monte Carlo simulations. It is thus worth some special treatment to reduce 
the computational load. In principle, the term R tJ appearing in the denominator of Eq. ([5]) 
should be strain-dependent. However, as we have chosen to expand the intersite interactions 
between local modes only up to harmonic order, it is consistent to ignore this effect, since 
strain-induced changes of the dipole-dipole interaction will enter only at higher order. This 
is equivalent to fixing the reciprocal lattice vectors G and all the atomic position vectors 
R,j. The dipole energy can then be written as 


Z 1 'y ( Qij,a/3'U j i,a'Uj,l3-i 
ij,a/3 


(7) 
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( 8 ) 


Here, a and f3 denote Cartesian components. The matrix Q is thus treated as a constant; 
it is calculated once and for all, and stored for later calculation of the dipole energy. This 
strategy increases the computational efficiency by at least one order of magnitude. 


D. Short-range interaction 

^short('| u |j th e energy contribution due to the short-range interactions between neigh¬ 
boring local modes, with dipole-dipole interaction excluded. This contribution stems from 
differences of short-range repulsion and electronic hybridization between two adjacent local 
modes and two isolated local modes. Together with the dipole-dipole interaction, this in¬ 
teraction determines the soft-mode energy away from the zone center. Expanded up to the 
second order, it can be written as, 

s*“({ U })dEEW^ P) 

Z i^j a/3 

The coupling matrix Ji h0 .d is a function of and should decay very fast with increasing 
Ry |. Here, we consider the short-range interaction up to third nearest neighbor (nn), 
whose local mode shares atoms with the local mode on the origin. Local modes between 
further neighbors involve displacements of atoms at least two hops away (in tight-binding 
language) and their core-core repulsion or hybridization should be much less important than 
the dipole-dipole interaction which is taken care of in £' dpl . 

The interaction matrix J,:y a /3 can be greatly simplified by symmetry. For a cubic lattice, 
we have, 

1st nn : J ijM3 = (j i + (j 2 - ji)\Rij, a \)S a p , 

2nd nn : J ijM3 = (j 4 + V2(j 3 ~ j 4 ,)\Rij, a \)$ a p 

- 0 * 2 .](1 $ a /3 ) j 

3rd 1111 . Jij,a@ j 6 *^a /3 T 8 / 7 /f;/. Rtj.3 (1 3 0 / 3 ) , (30) 
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where Rij >a is the a component of R ij/Rij. So we have only seven interaction parameters for 
a cubic lattice. The coefficients ji, j 2 , • jr in the above equations have physical meanings 
that are sketched schematically in Fig. |2j. For example, j\ represents the interaction strength 
of “ 7 r” like interactions between first neighbors. 

E. Elastic energy 

We will describe the state of elastic deformation of the BaTiOs crystal using local strain 
variables r]i (Rj), where the Voigt convention is used (/ = 1,..., 6) and Rj labels a cell center 
(Ti) site. In fact, the six variables per unit cell {^(R, ; )} are not independent, but are 
obtained from three independent displacement variables per unit cell. In our analysis, these 
are taken as the dimensionless displacements v(Rj) (in units of the lattice constant a) 
defined at the unit cell corner (Ba) positions Rj + (a/2, a/2, a/2). In terms of these, the 
inhomogeneous strain variables rjij (Rj) are defined in the next subsection. Because of our 
use of a periodic supercell in the Monte Carlo simulations, however, homogeneous strain 
deformations are not included in the configuration space {v(R ; )}. Therefore, we introduce 
six additional homogeneous strain components g^.i to allow the simulation cell to vary in 
shape. The total elastic energy is expanded to quadratic order as 

E elas ({m }) = Ef^ttrnA) + E^dmA), ( 11 ) 

where the homogeneous strain energy is simply given by 

E £ as ({vn,i}) = y5n(4 ;1 + 4,2 + 4,3) 

+ n B V2 (m,im,2 + + m,3Vn,i) 

N 

+ y £44 (4,4 + 4,5 + 4,6)' (12) 

Here B u, B u , and R 44 are the elastic constants expressed in energy units (Bn = a 3 Cn, 
etc.), and N is the number of primitive cells in the supercell. 

Rather than use an expression like (|T^) for the inhomogeneous strain energy, we have 
found it preferable to express this part directly in terms of the v(Rj).0 This approach keeps 
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the acoustic phonon frequencies well behaved throughout the Brillouin zone. To satisfy 
requirements of invariance under translations and rotations of the crystal as a whole, the 
energy is expanded in scalar products of differences between the v(Rj). The cubic crystal 
symmetry leads to a great reduction of the independent parameters in the expansion. The 
energies of the long-wavelength strain deformations can be reproduced by an expansion of 
the form 

E? las = ]T {711 MRi) - v x (Ri ± x )] 2 

i 

+ 712 MR.,;) - MR* ± X )]K(R-,) - VyiR, ± y)] 

T 744 \px (f^,) v x (bh ± y) + Uy(Rj) Vy (R, i x) ] 

+ cyclic permutations j (13) 

corresponding to bond stretching, bond correlation, and bond bending, respectively. Here, 
x = ax, y = ay, z = az, and ± indicates multiple terms to be summed. The 7 coefficients 
are related to the clastic constants by 711 = Bn/ 4, 712 = B 12 / 8 , and 744 = R 44 / 8 . 

F. Elastic-mode interaction 

To describe the coupling between the clastic deformations and the local modes, we use 
the on-site interaction 

M) = \ £ E Bmh(HiK(R>/s(Ri) ■ (u) 

2 i lot /3 

As a result of cubic symmetry, there are only three independent coupling constants Bi a p: 

B\xx E2yy B% zz , 

B\yy B ] zz B2xx B‘2zz B%XX B% yy , 

Bilyz B^zy B§XZ B 5 ZX B(jxy B(yyx ■ 

The strain contains both homogeneous and inhomogeneous parts, rji (Rj) = 7 H,z(R-i) + 
rji t i (Rj). The latter are expressed in terms of the local displacement vectors v as follows. 
We first define the six average differential displacements associated with site Rj as 
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Av xx = [w(Rj - d - x) - v x (R, - d)] , 

d=0,y,z,y+z 

A v xy = [v y (Ri - d - x) - i> y (Ri - d)] , 

d=0,y,z,y+z 

and their cyclic permutations, recalling that v(R;) is associated with position Rj + 
(a/2, a/2, a/2). Then 77^1 (R*) = Au xx / 4 , 771,4 (R*) = (Av yz + An 22/ )/ 4 , etc. 


III. FIRST-PRINCIPLES CALCULATIONS 


We have shown that, with the two approximations we made, the total energy functional 
of the perovskite system is fully specified by a set of parameters. These parameters can 
be obtained from first-principles calculations. We use density-functional theory within the 
local density approximation (LDA). The technical details and convergence tests of the cal¬ 
culations can be found in Refs. 0 ® The most important feature of the calculations is 
the Vanderbilt ultra-soft pseudopotential, 111 which allows a low energy cutoff to be used 
for first-row elements. This makes high-accuracy large-scale calculations of materials in¬ 
volving oxygen and 3d transition metal elements affordable. The ultra-soft pseudopotential 
also allows for exceptionally transferable pseudopotentials. It ensures that all-electron and 
pseudo-atom scattering properties agree over a very large energy range and preserves the 
chemical hardness of the atom. A generalized Kohn-Sham functional is directly minimized 
using a preconditioned conjugate gradient method.0011 We use a ( 6 , 6 , 6 ) Monkhorst-Pack 
k-point meshil for single-cell calculations, i.e., 216 k-points in the full Brillouin zone (FBZ). 
For supercell calculations, the k mesh is kept the same to minimize errors due to the k-point 
sampling. Therefore, a smaller number of k-points is used because of the smaller FBZ. 

We start with the construction of the local mode vectors. All the eigenvalues and eigen¬ 
vectors of the force-constant matrix at k = 0 for the cubic BaTiC >3 structure are calculated 
from frozen-phonon calculations, as in Ref. [U]. The mode with imaginary frequency is 
identified as the soft mode. The soft-mode eigenvector has been reported previously0 as 
^Ba = 0.20, £xi = 0.76, £o|| — —0.53, and £o_l = —0.21. The local mode is then constructed 
from it using the scheme proposed in II. 1 (Fig.|lj). 
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Determination of many of the parameters in the effective Hamiltonian involves only cal¬ 
culations for zone-center distortions. These parameters have been reported previously.01il 
They include the fourth-order terms of on-site energy a and 7 ; the elastic constants 
Bu, Bi 2 , -B 44 ; and the on-site elastic-mode interaction parameters Bi xx , Bi yy , B± yz . The 
mode effective charge Z* of Eq. (j4j) is calculated from the values Z* A = 2.75, Z^= 7.16 
Zq ,||=—5.69, and Z* ox =—2.11 published in Ref. [T^. (The resulting value Z*= 9.96 is slightly 
different from the one given in Ref. [12|; there, the eigenvector of the dynamical matrix, not 
the force-constant matrix, was used.) We use the experimental value = 5.24 of the op¬ 
tical dielectric constant, since for this quantity, the LDA seems not to be a well-justified 
approximation, while exact density-functional theory results are not accessible. We find, 
however, that the effect of a small inaccuracy in the dielectric constant affects thermody¬ 
namic properties such as transition temperatures only slightly. 

The second-order energy parameter k for zone-center distortions is a linear combination 
of the local mode self energy parameter the intersite interactions j t , and the dipole- 
dipole interaction. The calculations of intersite interaction parameters involve determination 
of the energy for distortions at the zone-boundary k-points X=( 7 r/a, 0 , 0 ), M=( 7 r/a, 7 r/a, 0 ), 
and R=( 7 r/a, 7 r/a, 7 r/a), where a is the lattice constant. Five frozen-phonon calculations 
on doubled unit cells are sufficient to extract all the information available from these k- 
points. The arrangement of the local mode vectors for each case, as well as for the zone- 
center distortion at T=(0,0,0), are shown in Fig. 0. The actual ionic configurations are 
constructed by superpositions of displacements associated with adjacent local modes. For 
example, letting u be the amplitude of the Ti-centered local mode defined in Sec. |I 1 A[ the 
displacement of the Ti atoms is just tdjri in Fig. 0(a) and itdyri in Fig. 0(b), while Ba is 
affected by eight neighboring local modes so that its displacement is 8 x uf Ba/8 = uf Ba in 
Fig. 0(a) and 0 in Fig. 0(b). 

The above five doubled-ccll calculations can be used to determine the parameters j 1 , 
j 2 , J 3 , J 4 , and je. The determination of j '5 + 2 j 7 requires a four-cell calculation involving 
20 atoms with low symmetry [Fig. 0(g)]. Table | lists the energy expressions for all the 
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configurations calculated in terms of the quadratic expansion parameters. 

A further decomposition of j 5 and j 7 would require an expensive eight-cell calculation. 
Furthermore, the interaction parameter j 7 is the third nearest-neighbor interaction and 
is thus presumably not very important. This argument is justified by our Monte Carlo 
simulations which show that the calculated transition temperature is insensitive to different 
decompositions of j '5 and j 7 . This prompts us to make an approximate decomposition 
based on a simple physical argument: we expect the interaction to be smallest for two 
interacting local modes oriented such that reversing the relative sign of the vectors produces 
the least change of bond lengths. Applied to third nearest neighbors, this argument implies 
j 6 — 2 j 7 = 0 , thus fixing the value of j 7 . 

The resulting interaction parameters are shown in Table p, together with other param¬ 
eters published previously. It may be surprising to see that the on-site k 2 is positive, while 
the cubic structure is known to be unstable against k — 0 distortion. The cubic structure 
is thus stable against forming an isolated local mode; instability actually comes from the 
intersite interactions between local modes. To be more precise, we find that it is the long- 
range Coulomb (dipole-dipole) interaction which makes the ferroelectric state favorable. If 
we turn off the dipole-dipole interaction by setting the effective charge Z* = 0, we find that 
the ferroelectric instability disappears. This is consistent with the previous point of view 
that long-range Coulomb forces favor the ferroelectric state, while short-range repulsions 
favor the nonpolar cubic state. 

From Table |TJ, we see that the intersite interaction parameters decay very fast with 
increasing distance, indicating the short-range nature of the intersite interactions after the 
long-range Coulomb interactions have been separated out. The ratio of the magnitudes of 
the strongest first-, second-, and third-neighbor interactions turns out to be approximately 
1 : 0.23 : 0.09. This decays even faster than the dipole-dipole interactions, for which the 
corresponding ratio (oc 1 /R 3 ) is 1 : 0.35 : 0.19. These results help justify our approximation 
of including only up to third nearest neighbors for the short-range interactions. 
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IV. MONTE CARLO SIMULATIONS 


For the quantitative study of non-universal finite-temperature behavior of a given model, 
Monte Carlo simulationElIl has emerged as the most reliable and powerful technique. It is 
especially appropriate for a model such as ours, with two continuous vector degrees of free¬ 
dom per unit cell and both short and long range interactions, for which analytical approaches 
snch as renormalization group or high-temperature expansions would be cumbersome and 
involve additional approximations. In comparison, the Monte Carlo approach requires only 
the ability to compute changes in total energy as the configuration is changed. Furthermore, 
with suitable analysis of statistical errors and finite-size effects, the results of Monte Carlo 
simulation can be made arbitrarily accurate. Finally, with little additional effort, a number 
of physical quantities can be computed to aid in characterization of the transition. 

We solve the effective Hamiltonian [Eqs. (H), (|7|), (0h (HHb and ([14])] using Monte Carlo 
simulations with the Metropolis algoritlimil on an L x L x L cubic lattice with periodic 
boundary conditions. Since most energy contributions (except E dpl ) are local, we choose 
the single-flip algorithm. That is, a trial move consists of an attempted update of a single 
variable, after which the total energy change is calculated to determine whether to accept 
the move. The step sizes are adjusted to ensure an acceptance ratio of approximately 0.2. In 
one Monte Carlo sweep (MCS), we first make a trial move on each Uj in sequence, then each 
Vj in sequence, then iterate several times (typically 2 L times) on the homogeneous strain 
variables. For L = 12, each MCS takes about one second on an HP 735 workstation. The 
typical correlation time for the total energy is found to be several hundred MCS close to the 
phase transition; this long correlation time makes certain new MC techniques using energy 
distribution functional unfavorable. The correlation times for the local mode amplitudes 
are one order of magnitude shorter, and thus 10,000 MCS are usually enough to equilibrate 
and to obtain averages of local mode variables with statistical error < 10%. 

In our simulation, we concentrate on identifying the succession of low-temperature 
phases, determining the phase transition temperatures, and extracting qualitative features 
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of the transitions. This analysis will allow us to identify the features of the effective Hamil¬ 
tonian which most strongly affect the transition properties. For these purposes, it is most 
convenient to monitor directly the behavior of the homogeneous strain and the vector order 
parameter. In the case of the ferroelectric phase transition, the latter is just the average 
local-mode vector u = Uj/iV, which is proportional to the polarization. Here, u* is the lo¬ 
cal mode vector at site i and N is the total number of sites. As a reference, the average local 
mode amplitude u = |uj|/iV is also monitored. To avoid effects of symmetry-equivalent 
rotations of the order parameter and to identify the different phases clearly, we accumulate 
the absolute values of the largest, middle, and smallest components of the averaged local¬ 
mode vector for each step, denoted by rq, w 2 , and u 3 , respectively (tq > u 2 > u 3 ). The cubic 
(C), tetragonal (T), orthorhombic (O), and rhombohedral (R) phases are then characterized 
by zero, one, two, and three non-zero order-parameter components, respectively. The effect 
of symmetry-equivalent rotations on the homogeneous strain is handled analogously, with 
the largest, the medium and the smallest linear strain components denoted by 771, 772 and 773, 
respectively, and shear strain components by 774, 775 and 7 ] e . 

The ferroelectric phase transition is very sensitive to hydrostatic pressure, or equivalently, 
to lattice constant. The LDA-calculated lattice constants are typically 1% too small, and 
even this small error can lead to large errors in the zero-pressure transition temperatures. 
One approach which largely compensates for the effect of this systematic error is to exert a 
negative pressure that expands the lattice constant to the experimental value. We determine 
the value of the pressure by calculating volumes for four different phases and comparing 
with experimental measurements. 0 We find P = —4.8 GPa gives the best overall agreement 
(although the application of pressure does lead to a slight change in the low temperature 
structure). Except for the simulations for the construction of the temperature-pressure phase 
diagram, the following simulations and analysis are for this pressure. 
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V. RESULTS AND DISCUSSION 


In this section, the finite-temperature behavior of the model is presented and analyzed. 
First, we examine the order parameters as a function of temperature in a typical simulation 
to obtain a measure of the transition temperatures. From the results of simulations for a 
range of pressures, we construct the temperature-pressure phase diagram. For the system 
at ambient pressure, more detailed simulations are performed. The order of the transitions, 
the nature of the paraelectric phase and the properties of the low-temperature phases are 
investigated and compared with experimental observations. Finally, the roles played by 
different terms in the effective Hamiltonian and the sensitivity of the results to various 
approximations are examined. 

A. Order parameters and phase diagram 

We start the simulations at a high temperature (T > 400K) and equilibrate for 10,000 
MCS. The temperature is then reduced in small steps, typically 10K. After each step, the 
system is allowed to equilibrate for 10,000 MCS. The order parameter averages are then accu¬ 
mulated over a period of 7,000 MCS, yielding a typical standard deviation of less than 10%. 
The temperature step size is reduced and the number of MCS is doubled for temperatures 
close to the phase transition. 

Fig. [| shows the averaged local-mode vectors iq, w 2 , u 3 and averaged local mode ampli¬ 
tude u as functions of temperature in a typical simulation for an L = 12 lattice at P = —4.8 
GPa. At high temperatures, u 3 , w 2 , and u 3 are all very close to zero. As the system is 
cooled down below 295K, tq increases and becomes significantly larger than tt 2 or u 3 . This 
indicates the transition to the tetragonal phase. Two additional phase transitions occur 
as the temperature is reduced further. The T-0 transition (sudden increase of w 2 ) occurs 
at 230K and the O-R transition (sudden increase of u 3 ) occurs at 190K. The sequence of 
transitions exhibited by the simulation is the same as that observed experimentally. 
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The averaged homogeneous strain variables obtained from the above simulation are 
shown in Fig. [|. These strains are measured relative to the LDA calculated equilibrium 
cubic structure, so the linear strains are significantly non-zero at higher temperatures due 
to the negative pressure applied. As expected, the simulation cell changes shape at the 
same temperatures at which the jumps of order parameter components are observed. At 
high temperatures, we have approximately rj i = 772 = 77 3 and 774 = 775 = ? 7 6 = 0, correspond¬ 
ing to the cubic structure. As the system is cooled down, the shape of the simulation cell 
changes to T, O, and R phases. The orthorhombic (O) structure has a non-zero shear strain, 
in agreement with the centered orthorhombic structure observed experimentally. Quantita¬ 
tively, our calculated distortions are also in good agreement with the experiment, with the 
calculated distortions tending to be slightly smaller. For example, 774 — 773 for the tetragonal 
phase is 1.1% as measured from experimental and 0.9% from our calculation. 

The simulations are repeated for a range of applied pressures to obtain the temperatures 
at which the order parameter components and homogeneous strain jump on cooling down. 
The resulting temperature-pressure phase diagram is shown in Fig. |S]. (This measure of 
the transition temperature is actually a lower bound, due to hysteresis around 5% for T-0 
and O-R transitions and negligible for the C-T transition, to be discussed further below.) 
All three transition temperatures decrease almost linearly with increasing pressure. At the 
experimental lattice constant, the values for dTjdP are found to be —28, —22, and —15 
K/GPa for the C-T, T-O, and O-R transitions respectively. The experimental values for 
the C-T transition range from —40 K/GPail to —66 K/GPa.i! For the T-0 transition the 
measured value is —28 K/GPa^l and for the O-R transition it is —10 to —15 K/GPa.0 At 
pressures as high as P = 5GPa, the sequence of phases C-T-O-R is still observed in the 
simulation. When the pressure is increased further, the phase boundary in the simulation 
becomes unclear due to fluctuations. Our calculated critical pressure (beyond which the 
cubic structure is stable at T=0 K) is P c = 8.4GPa. Taking into account the pressure 
correction for the LDA volume underestimate, this corresponds to a predicted physical 
P c = 13.2 GPa. We are not aware of any experimental value for P c with which to compare 
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this value. However, we find that the magnitude of our dT c /dP is significantly smaller than 
experimental value, at least for the C-T and T-0 transitions. This may partly be due to the 
neglect of higher-order strain coupling terms in the effective Hamiltonian. We have tested 
the effect of including a volume dependence for the short-range interaction parameters j t . 
This correction does not change the sequence of phases, and it only increases the magnitudes 
of dT c jdP slightly. Therefore, our results are reported without this correction. The accuracy 
of the phase diagram may be further improved by including higher-order terms in the elastic 
energy, or the coupling of j t to anisotropic strain. 


B. Hysteresis and latent heat 


For the investigation of the order of the transitions, the nature of the paraelectric phase 
and the properties of the low-temperature phases, we performed more detailed simulations 
at P = —4.8GPa for the three system sizes L = 10, 12 and 14. In the cooling-down 
simulations, the length of each simulation was increased from 10,000 to up to 35,000 MCS 
at temperatures close (±10K) to the phase transition to include a longer equilibration. The 
size of the temperature step was decreased to 5K or less in the vicinity of the transition. In 
addition, a heating simulation was performed, starting from the lower-temperature phase, 
to detect any possible hysteresis. The calculated transition temperatures, obtained as the 


average of the cooling and heating transition temperatures, are given in Table |TTT| . The 
error estimates in the table are determined by the width of the hysteresis, which persists 
even for the longest simulation lengths considered. (The C-T transition temperature for 
L = 10 is difficult to identify because of large fluctuations between phases.) The calculated 
transition temperatures are well converged with respect to system size, and are in reasonable 
agreement with experiment. Table 0 also gives the saturated spontaneous polarization p s 
at T = 0 in the R phase, and just above the O-R and T-0 transitions in the O and T phases 
respectively. These are calculated from the average local mode vector u and the local mode 
Z*. We find almost no finite-size effect for this quantity, and the agreement with experiment 
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is very good for the O and T phases. The disagreement for the R phase may result in part 
from twinning effects in the experimental sample.!*] 

From the jumps in structural parameters and the observed hysteresis in heating and 
cooling, we conclude that the phase transitions are first order. An accurate determination 
of the latent heats would require considerable effort;ifl here, we only try to provide estimates 
sufficiently accurate for meaningful comparison with experiment. We approach each transi¬ 
tion from both high-temperature and low-temperature sides until the point is reached where 
both phases appear equally stable. (That is, the typical time for the system to fluctuate into 
the opposite phase is roughly independent of which phase the simulation is started in.) The 
difference of the average total energy is then the latent heat.il This approach is practical 
as long as some hysteresis is present. The calculated latent heats (Table B show non- 
negligible hnite-size dependence. Taking this into account, we find that the latent heats for 
all three transitions are significantly non-zero and in rough agreement with the rather scat¬ 
tered experimental data. For the T-0 and O-R transitions, the first-order character of the 
transition is predicted by Landau theory, since in these two cases the symmetry group of the 
low temperature structure is not a subgroup of that of the high temperature structure. For 
the C-T transition, the first-order character is not a consequence of symmetry, but rather of 
the values of the effective Hamiltonian parameters. Although it has the largest latent heat 
of the three transitions, the C-T transition also exhibits large hnite-size effects in the latent 
heat and in the smearing of order parameter components and strain discontinuities in the 
simulation (Figs. [|,|5]). This suggests the presence of long-wavelength fluctuations character¬ 
istic of second-order phase transitions and critical phenomena, and the classification of the 
C-T transition as a weak first-order transition. 

C. Displacive vs. order-disorder 

Using the microscopic information available in the simulations, we are able to investigate 
the extent to which the cubic-tetragonal transition can be characterized as order-disorder 
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or displacive. These possibilities can be distinguished by inspecting the distribution of the 
real-space local-mode vector Uj in the cubic phase just above the transition. A displacive 
(microscopically nonpolar) or order-disorder (microscopically polar) transition should be 
characterized by a single-peaked or double-peaked structure, respectively. The distribution 
of u x at T = 400K is shown in Fig. 0. It exhibits a rather weak tendency to a double-peaked 
structure, indicating a transition which has some degree of order-disorder character. We also 
see indications of this in the u—T relation in Fig. f|. Even in the cubic phase, the average of 
the local-mode magnitude u is significantly non-zero and close to that of the rhombohedral 
phase, while the magnitudes of the average local mode components change dramatically 
during the phase transitions. 

In reciprocal space, a system close to a displacive transition should show large and 
strongly temperature-dependent fluctuations of certain modes associated with a small por¬ 
tion of the Brillouin zone (BZ) (for a ferroelectric transition, near T). For an extreme order- 
disorder transition, on the other hand, the fluctuations are expected to be distributed over 
the whole BZ. For BaTi0 3 , we calculated the average Fourier modulus F(k, T) = (|w(k)| 2 ) 
for eigenmodes at several high-symmetry k-points (along T-X, T-M and T-R) for a range 
of temperatures above the C-T transition. These eigenmodes are identified by their sym¬ 
metry properties as one longitudinal optical (LO) branch and two transverse optical (TO) 
branches at each point. For a purely harmonic system, T/F(k,T) can be shown to be a 
temperature-independent constant proportional to the square of the eigenfrequency cn 2 (k) 
of the corresponding eigenmode. A strong decrease of T/F(k,T ) as T —> T c from high 
temperature can be interpreted as mode softening due to anharmonicity. 

The results at the k-point (7 t/4 a, ir/Aa, 0) illustrate the main features of this analysis. 
From symmetry, three nondegenerate eigenmodes LO, TOl, T02 are identified. The polar¬ 
ization of LO, TOl, and T02 are in the direction of x + y, x — y, and z, respectively. For 
each eigenmode, the temperature dependence of the calculated u; 2 (k, T) is shown in Fig. | 8 ]. 
The almost linear behavior of cn 2 (k, T) vs. T (the Curie-Weiss form) is observed for the other 
k-points as well. Both the LO and TOl branches are almost temperature independent. The 
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T02 branch is strongly temperature dependent and is thus a “soft” mode. According to the 
soft-mode theory of structural phase transitions, T c is the lowest temperature at which all 
ta 2 (k, T) > 0. Linear extrapolation indicates that the T02 mode frequency goes to zero at 
T 200K, which is a lower bound for T c , consistent with the value obtained in Monte Carlo 
simulations. A similar calculation of cn 2 (k, T) for the TO modes at T=(0,0,0) extrapolates 
to zero at the higher temperature T « 300K, in excellent agreement with the Monte Carlo 
value of 295K. 

Within this formalism, the microscopic character of the paraelectric phase is determined 
by the extent of the soft mode in the BZ. We define a quantity 


2F(k, 350K) 
^ ’ ~ F{ k, 700K) 


(15) 


to indicate the hardness of the modes. In Fig. |9], p(k) is shown for the various k-points 
along some special directions in the BZ. If p(k) < 1, the corresponding eigenfrequency 
extrapolates to zero at some positive temperature, and the mode is regarded as soft. If 
u; 2 (k) is independent of temperature, p(k) = 2, corresponding to the hardest mode. 

For all the k-points considered, all the LO modes are almost temperature independent 
[p(k) = 2] and are not included in the figure. Along the T-X direction, the doubly-degenerate 
TO modes are soft all the way to the zone boundary. In contrast, along the T-R, direction, 
both TO modes become hard immediately after leaving the T point. Along the T-M di¬ 
rection, the TOl mode becomes hard quickly, while the T02 branch remains soft at least 
halfway to the zone boundary. This behavior, especially along T-X, does not conform com¬ 
pletely to the displaeive limit. This supports the interpretation of the C-T transition as 
intermediate between displaeive and order-disorder, with a slightly stronger order-disorder 
character. Thus, from the example of BaTi0 3 , it seems that a positive on-site quadratic 
coefficient does not automatically imply a displaeive character for the transition. Rather, 
the relevant criterion is the extent to which the unstable modes extend throughout the BZ. 


23 



D. Roles of different interactions 


Our theoretical approach allows us to investigate the roles played by different types of 
interaction in the phase transition. First, we study the effect of strain. Recall that the 
strain degrees of freedom were separated into local and homogeneous parts, representing 
finite- and infinite-wavelength acoustic modes, respectively. Both parts were included in the 
simulations. If we eliminate the local strain (while still allowing homogeneous strain), we 
find almost no change in the transition temperatures. This indicates that the effect of the 
short-wavelength acoustic modes may not be important for the ferroelectric phase transition. 
If the homogeneous strain is frozen at zero, however, we find a direct cubic-rhombohedral 
phase transition, instead of the correct series of three transitions. This demonstrates the 
important role of homogeneous strain. 

Second, we studied the significance of the long-range Coulomb interaction in the simu¬ 
lation. To do this, we changed the effective charge of the local mode (and thus the dipole- 
dipole interaction strength), while modifying the second-order self-energy parameter k-z , so 
that the frequencies of the zone-center and zone-boundary modes remain in agreement with 
the LDA values. We found only a slight change (10%) of the transition temperatures when 
the dipole-dipole interaction strength was doubled. However, elimination of the dipole- 
dipole interaction altogether changed the results dramatically; the ground state becomes a 
complex antiferroclectric structure similar to the room-temperature structure of PbZrOs. 
This result shows that it is essential to include the long-range interaction, although small 
inaccuracies in the calculated values of the effective charges or dielectric constants may not 
be very critical. 

Third, we investigated the sensitivity of our results to variations of the short-range inter¬ 
action parameters. We find the accuracy of the first-neighbor interaction parameters (jij' 2 ) 
is very important, and a mere 10% deviation can change the calculated transition temper¬ 
atures dramatically, and can sometimes even change the ground state structure. Second 
nearest-neighbor interactions are less important, and for the third-neighbor interactions, 
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even a 100% change does not seem to have a strong effect on the values of T c . This result 
is to be expected, and partly justifies our choice of including only up to third neighbors for 
the short range interactions. We have also tested the effect of our assumption j 6 — 2 j 7 = 0. 
We find that any reasonable choice leads to a barely noticeable change in T c . 

In short, highly accurate LDA calculations do appear to be a prerequisite for an accurate 
determination of the transition temperatures, but as long as certain features of the energy 
surface are correctly described, other approximations can be made without significantly 
affecting the results. 


VI. CONCLUSIONS 

We have developed a first-principles approach to the study of structural phase transi¬ 
tions and finite-temperature properties in perovskite compounds. We construct an effective 
Hamiltonian based on Taylor expansion of the energy surface around the cubic structure, 
including soft optical modes and strain components as the possible distortions. The ex¬ 
pansion parameters are determined by first-principles density-functional calculations using 
Vanderbilt’s ultra-soft pseudopotential. 

We have applied this scheme to BaTi0 3 and calculated the pressure-temperature phase 
diagram. We have obtained the sequence of low-temperature phases, the transition temper¬ 
atures, and the spontaneous polarizations, and found them to be in good agreement with 
experiment. We find that long-wavelength acoustic modes and long-range dipolar interac¬ 
tions both play an important role in the phase transition, while short-wavelength acoustic 
modes are not as significant. Accurate LDA calculations are required for accurate deter¬ 
mination of the transition temperatures. The C-T phase transition is not found to be 
well described as a simple displacive transition; on the contrary, if anything it has more 
order-disorder character. 

With slight modifications, our approach should be applicable to other perovskite com¬ 
pounds, as long as their structure is close to cubic and there are some low energy distortions 


25 



responsible for the phase transitions. It can be easily applied to ferroelectric materials like 
PbTi0 3 (under study by another groupEI) and KNb0 3 . It can also be applied to antiferro- 
electric materials like PbZr0 3 . The application to antiferrodistortive materials like SrTi0 3 
is slightly more involved, though also successful.il 
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TABLES 

TABLE I. The energy per 5-atom unit cell (excluding the dipole energy) in terms of intersite 
interaction parameters of Fig. 2, for configurations shown in Fig. 3. 


Configuration 

expression 

(a) 

+ 2 ji + j 2 + 4j 3 + 2j 4 + 4j 6 

(b) 

K2 + 2,71 - ,72 - 4j 3 + 2,74 - 4jg 

(c) 

«2 + .72 - 2J4 - 4.76 

(d) 

«2 - 2ji + ,72 - 4.73 + 2j 4 + 4j 6 

(e) 

K2 + .72 - 2,74 + 4j 6 

(f) 

«2 - 2j'i - ,72 + 4.73 + 2j 4 - 4j 6 

(g) 

K2/2 + Ji - 2,75 - 4,77 


TABLE II. Expansion parameters of the Hamiltonian for BaTiOs- Energies are in Hartrees. 


On-site 

«2 

0.0568 

a 

0.320 

7 

-0.473 


h 

-0.02734 

32 

0.04020 



Intersite 

J3 

0.00927 

34 

-0.00815 

k 

0.00580 


k 

0.00370 

37 

0.00185 



Elastic 

B u 

4.64 

B\2 

1.65 

-R 44 

1.85 

Coupling 

B\xx 

-2.18 

B\yy 

-0.20 

B/±yz 

-0.08 

Dipole 

z* 

9.956 

^OO 

5.24 
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TABLE III. Calculated transition temperatures T c , saturated spontaneous polarizations p s 


and estimated latent heats l, as a function of simulation cell size. 



phase 

O 

T—1 

II 

L = 12 

L = 14 

expt a 


O-R 

197±3 

200 T10 

200±5 

183 

T c (K) 

T-0 

230T10 

232±2 

230T10 

278 


C-T 

~290 

296±1 

297±1 

403 


R 

0.43 

0.43 

0.43 

0.33 

p s (C/m 2 ) 

0 

0.35 

0.35 

0.35 

0.36 


T 

0.28 

0.28 

0.28 

0.27 


O-R 

50 

60 

60 

33-60 

l (J/mol) 

T-0 

90 

90 

100 

65-92 


C-T 

- 

- 

150 

196-209 


a Ref. H 
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FIGURES 

FIG. 1. The structure of the cubic perovskite compound BaTiC^. Ba, Ti and O atoms are 
represented by shaded, solid, and empty circles respectively. Lengths of the vectors indicate the 
relative magnitudes of the displacements for a local mode polarized along x. 

FIG. 2. The independent intersite interactions corresponding to the parameters j\ , U (first 
neighbor), j' 3 , j' 4 , j '5 (second neighbor), and jg and j 7 (third neighbor). 

FIG. 3. The local mode arrangements for which first-principles total-energy calculations were 
performed to determine the intersite interaction parameters. The arrangements can be labeled by 
the wavevector k and a polarization vector (p). The arrows represent local mode vectors. The 
dotted lines indicate the unit cells of the simple cubic lattice. The solid lines show the supercells 
used in the calculations, (a) T, p = z ; (b) X, p = z ; (c) X, p = x ; (d) M, p = z ; (e) M, p = x 
; (f) R, p = z ; (g) four-cell calculation. 

FIG. 4. The averaged largest, middle, and smallest components u\, 112 , 113 and amplitude u of 
local modes as a function of temperature in the cooling-down simulation of a 12 x 12 x 12 lattice 
described in Section IV. The dotted lines are guides to the eye. 

FIG. 5. The averaged homogeneous strain i]h as a function of temperature in the cooling-down 
simulation of a 12 x 12 x 12 lattice described in Section IV. The strains are measured relative to the 
LDA minimum-energy cubic structure with lattice constant 7.46 au. The dotted lines are guides 
to the eye. 

FIG. 6 . The calculated pressure-temperature phase diagram. The cubic-tetragonal (C-T), 
tetragonal-orthorhombic (T-O) and orthorhombic-rhombohedral (O-R) transitions are labeled by 
solid triangles, circles and squares, respectively. The vertical dash-dot line at P=— 4.8 GPa corre¬ 
sponds to zero pressure in the experiment to compensate for the LDA volume error. 

FIG. 7. The probability distribution of the Cartesian component of the local mode variable u x 
in the cubic phase at T = 320K (solid line), 350K (dashed line), and 500K (dotted line). 
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FIG. 8. Temperature dependence of squared eigenfrequency uj 2 at k = (7r/4a, 7r/4a, 0) for (a) 


LO, (b) TOl, and (c) T02 modes. 


FIG. 9. Calculated mode hardness quantity p(k), Eq. (0), along special directions in the 
Brillouin zone. 
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